library(here)
setwd(here())

load(file = "data/intermediate/models_multivariatefmm.Rdata")
load("data/intermediate/pf.rda")
source("scripts/promise_functions.R")
library(mellonMisc)
library(brms)
library(ggplot2)
promise.index$text[!promise.index$text %in% pf$Snippets]
promise.index$text[promise.index$text=="Keep the Trident nuclear deterrant" ] <-  "Keep the Trident nuclear deterrent" 
promise.index$text[promise.index$text== "Create a 'Veterans Board' to co-ordinate treatment of military veterans" ] <-   "Create a 'Veterans Board' to coordinate treatment of military veterans"

pf <- pf[match(promise.index$text, pf$Snippets), ]

##### Fixed effects results for comparison
library(readxl)
library(haven)
library(dplyr)
fe <- read_stata("data/intermediate/fixed effects model weights.dta")
fe <- fe %>% slct(promise, importFixedAll = importAll, 
                  approvalFixedAll = approvalAll, 
                  importFixedCon = importCon, approvalFixedCon= approvalCon,
                  fixedSD) 

fe$importFixedAll <- fe$importFixedAll - mean(fe$importFixedAll)
fe$importFixedCon <- fe$importFixedCon - mean(fe$importFixedCon)

comb.sd.fe <- mean(c(sd(fe$importFixedAll), sd(fe$importFixedCon)))

fe$importWtAll <- prop.table(exp(fe$importFixedAll /  comb.sd.fe))
fe$importWtCon <- prop.table(exp(fe$importFixedCon /  comb.sd.fe))
fe <- safemerge(fe, pf, by.x = "promise", by.y = "id")


prop.table(tapply(fe$importWtAll, list(sq = fe$status_quo, status = fe$statuscomb), sum), 1)

prop.table(tapply(fe$importWtCon, list(sq = fe$status_quo, status = fe$statuscomb), sum), 1)





# re.all <- ranef(fit.all)
sd.samples <- posterior_samples(fit.all.fmm, "sd_")
hist(sqrt(sd.samples$sd_mmg1g2g3g4__import_Intercept ^ 2 + 
            sd.samples$sd_mmmii1mii2mii3mii4__import_Intercept ^ 2))

# plot(fit.all)
samples <- posterior_samples(fit.all.fmm, "mmg1g2g3g4|mmmii1mii2mii3mii4")
imp.samples <- samples[, paste0("r_", "mmg1g2g3g4__import", "[", rownames(promise.index), ",Intercept]")]
app.samples <- samples[, paste0("r_", "mmg1g2g3g4__approval", "[", rownames(promise.index), ",Intercept]")]

mii.imp.samples <- samples[, paste0("r_", "mmmii1mii2mii3mii4__import", "[", mii.lookup$code, ",Intercept]")] 
mii.app.samples <- samples[, paste0("r_", "mmmii1mii2mii3mii4__approval", "[", mii.lookup$code, ",Intercept]")] 

promise.index$mii <- pf$MII_cat1[match(promise.index$text, pf$Snippets)]
promise.index$mii.code <- mii.lookup$code[match(promise.index$mii, mii.lookup$cat)]

imp.samples <- imp.samples + mii.imp.samples[, promise.index$mii.code] 

class.samples <- posterior_samples(fit.all.fmm, "classprop")

fe.add <- t(sapply(class.samples$b_import_classprop,
                   function(x) x * group_assign$majorPromiseProb))

imp.samples <- imp.samples + fe.add

imp.samples <- imp.samples - matrix(inverse.rle(list(values = rowMeans(imp.samples), 
                                                     lengths = rep(dim(imp.samples)[2], 
                                                                   dim(imp.samples)[1]))), 
                                    byrow = T, ncol = 257)
library(matrixStats)

app.samples <- app.samples + mii.app.samples[, promise.index$mii.code] 
library(matrixStats)
library(mellonMisc)

rowsds <- rowSds(as.matrix(imp.samples))
mean(rowsds)

med.eff <- colMedians(as.matrix(imp.samples))

samples.c <- posterior_samples(fit.con.fmm, "mmg1g2g3g4|mmmii1mii2mii3mii4")

imp.samples.c <- samples.c[, paste0("r_", "mmg1g2g3g4__import", "[", rownames(promise.index), ",Intercept]")]
imp.samples.c.vc <- samples.c[, paste0("r_", "mmg1g2g3g4__import", "[", rownames(promise.index), ",votecon]")]
imp.samples.c <- imp.samples.c + imp.samples.c.vc

mii.imp.samples.c <- samples.c[, paste0("r_", "mmmii1mii2mii3mii4__import", "[", mii.lookup$code, ",Intercept]")]
mii.imp.samples.c.vc <- samples.c[, paste0("r_", "mmmii1mii2mii3mii4__import", "[", mii.lookup$code, ",votecon]")]
mii.imp.samples.c <- mii.imp.samples.c + mii.imp.samples.c.vc


app.samples.c <- samples.c[, paste0("r_", "mmg1g2g3g4__approval", "[", rownames(promise.index), ",Intercept]")]
app.samples.c.vc <- samples.c[, paste0("r_", "mmg1g2g3g4__approval", "[", rownames(promise.index), ",votecon]")]

app.samples.c <- app.samples.c + app.samples.c.vc

mii.app.samples.c <- samples.c[, paste0("r_", "mmmii1mii2mii3mii4__approval", "[", mii.lookup$code, ",Intercept]")] 
mii.app.samples.c.vc <- samples.c[, paste0("r_", "mmmii1mii2mii3mii4__approval", "[", mii.lookup$code, ",votecon]")] 

mii.app.samples.c <- mii.app.samples.c + mii.app.samples.c.vc

app.samples.c <- app.samples.c + mii.app.samples.c[, promise.index$mii.code] 
imp.samples.c <- imp.samples.c + mii.imp.samples.c[, promise.index$mii.code] 

class.samples.c <- posterior_samples(fit.con.fmm, "classprop")

fe.add.c <- t(sapply(class.samples.c$b_import_classprop,
                     function(x) x * group_assign$majorPromiseProb)) + 
  t(sapply(class.samples.c$`b_import_votecon:classprop`,
           function(x) x * group_assign$majorPromiseProb)) + 
  t(sapply(class.samples.c$b_import_classprop.c,
           function(x) x * group_assign$majorPromiseProb.c)) + 
  t(sapply(class.samples.c$`b_import_votecon:classprop.c`,
           function(x) x * group_assign$majorPromiseProb.c))


imp.samples.c <- imp.samples.c + fe.add

imp.samples.c <- imp.samples.c - matrix(inverse.rle(list(values = rowMeans(imp.samples.c), 
                                                         lengths = rep(dim(imp.samples.c)[2],
                                                                       dim(imp.samples.c)[1]))), 
                                        byrow = T, ncol = dim(imp.samples.c)[2])

rowsds.c <- rowSds(as.matrix(imp.samples.c))
comb.sd <- mean(c(rowsds, rowsds.c))

save(comb.sd, file = "data/intermediate/comb_sd.rda")

imp.z <- imp.samples / comb.sd

imp.exps <- apply(imp.z, MARGIN = 1, exp)
imp.exps<- t(imp.exps)
imp.props <- prop.table(as.matrix(imp.exps), 1)

imp.c.z <- imp.samples.c / comb.sd

imp.exps.c <- apply(imp.c.z, MARGIN = 1, exp)
imp.exps.c <- t(imp.exps.c)
imp.props.c <- prop.table(as.matrix(imp.exps.c), 1)



mii.sums <- miiSumCalc(imp.props)
mii.sums.c <- miiSumCalc(imp.props.c)

mii.comb.vals <- dtf(Issue = names(mii.sums), 
                     All = colMedians(as.matrix(mii.sums)), 
                     Con = colMedians(as.matrix(mii.sums.c)))
mii.comb.vals$All <- prop.table(mii.comb.vals$All)
mii.comb.vals$Con <- prop.table(mii.comb.vals$Con)

unwt.iss <- prop.table(table(pf$MII_cat1))
iss.name.lu <- dtf(dirty = sort(mii.comb.vals$Issue), clean = names(table(pf$MII_cat1)))
mii.comb.vals$Issue <- iss.name.lu$clean[match(mii.comb.vals$Issue, iss.name.lu$dirty)]
mii.comb.vals$Unweighted <- unwt.iss[mii.comb.vals$Issue]

write.csv(mii.comb.vals, file = "tables/mii.comb.vals.csv", row.names = F)

mii.plot.noeu <- prepMIIPlot(mii.sums, iss.exclude = "europe")
mii.plot.c.noeu <- prepMIIPlot(mii.sums.c, iss.exclude = "europe")

mii.plot <- prepMIIPlot(mii.sums)
mii.plot.c <- prepMIIPlot(mii.sums.c)
library(ggplot2)

joy.all.iss <- ggplot(data = mii.plot, aes(x = Proportion, y = Issue,
                                           fill = factor(fill.count)) ) + 
  scale_fill_manual(values = c("#3aa1fc", "#0b4e89")) +
  geom_density_ridges(rel_min_height = 0.01, 
                      scale = 4,   alpha = 1, colour = "white") + 
  scale_x_continuous(labels = percent, 
                     limits = c(0, max(c(mii.plot.c$Proportion, mii.plot$Proportion)))) + 
  xlab("% of Manifesto (posterior density)") + 
  ylab("") + theme_minimal() + guides(fill = F, colour = F) +
  theme(panel.grid.major.y = element_line(size=0.4, 
                                          color = "gray", linetype = 3),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.ticks.x = element_line(size=0.1, color = "black")) + 
  geom_segment(aes(x=Unweighted, xend = Median, y = Issue, yend = Issue, colour = Diff),
               arrow = arrow(length = unit(0.2, "cm")), size = 0.5) + 
  geom_point(aes(x = Unweighted, y = Issue), colour = "black", alpha = 0.1)  + 
  scale_color_gradient(low = "#dd250d", high = "#18a306")


joy.all.iss.noeu <- ggplot(data = mii.plot.noeu, aes(x = Proportion, y = Issue,
                                                     fill = factor(fill.count)) ) + 
  scale_fill_manual(values = c("#3aa1fc", "#0b4e89")) +
  geom_density_ridges(rel_min_height = 0.01, 
                      scale = 4,   alpha = 1, colour = "white") + 
  scale_x_continuous(labels = percent, 
                     limits = c(0, max(c(mii.plot.c$Proportion, mii.plot$Proportion)))) + 
  xlab("% of Manifesto (posterior density)") + 
  ylab("") + theme_minimal() + guides(fill = F, colour = F) +
  theme(panel.grid.major.y = element_line(size=0.4, 
                                          color = "gray", linetype = 3),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.ticks.x = element_line(size=0.1, color = "black")) + 
  geom_segment(aes(x=Unweighted, xend = Median, y = Issue, yend = Issue, colour = Diff),
               arrow = arrow(length = unit(0.2, "cm")), size = 0.5) + 
  geom_point(aes(x = Unweighted, y = Issue), colour = "black", alpha = 0.1)  + 
  scale_color_gradient(low = "#dd250d", high = "#18a306")

joy.all.iss.c <- ggplot(data = mii.plot.c, aes(x = Proportion, y = Issue,
                                               fill = factor(fill.count))) + 
  scale_fill_manual(values = c("#3aa1fc", "#0b4e89")) +
  geom_density_ridges(rel_min_height = 0.01, 
                      scale = 4,  colour = "white", alpha = 1) + 
  scale_x_continuous(labels = percent, 
                     limits = c(0, max(c(mii.plot.c$Proportion, mii.plot$Proportion))) ) + 
  xlab("% of Manifesto (posterior density)") + 
  ylab("") + theme_minimal() + guides(fill = F, colour = F) +
  theme(panel.grid.major.y = element_line(size=0.4, 
                                          color = "gray", linetype = 3),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.ticks.x = element_line(size=0.1, color = "black")) + 
  geom_segment(aes(x=Unweighted, xend = Median, y = Issue, yend = Issue, colour = Diff),
               arrow = arrow(length = unit(0.2, "cm")), size = 0.5) + 
  geom_point(aes(x = Unweighted, y = Issue), colour = "black", alpha = 0.1)  + 
  scale_color_gradient(low = "#dd250d", high = "#18a306")
library(gridExtra)
joy2 <- arrangeGrob(joy.all.iss + ggtitle("All respondents"), 
                    joy.all.iss.c + ggtitle("2017 Conservatives"), nrow = 2)
plot(joy2)

pf$cat_small <- pf$agreed_cat_string
pf$cat_small[!pf$cat_small %in% c("europe", "immigration", "defence", "health", 
                                  "education", "austerity", "pensions/ageing")] <- "other"

pf$Done <- as.numeric(pf$statuscomb=="Fulfilled")

prop.table(table(pf$statuscomb, pf$status_quo), 2)
median(rowSums(propDone(imp.props = imp.props, "no", pf = pf)))

props.total <- dtf(Complete = rowSums(propDone(imp.props, sq = c("yes", "no"), pf= pf)))

unweighted <- prop.table(table(pf$Done))["1"]
density.all <- ggplot(data = props.total, aes(x = Complete)) + 
  # scale_fill_manual(values = c("#3aa1fc", "#f49241")) +
  geom_density(colour = "#568aa3", alpha = 0.4, fill = "#3aa1fc") + 
  scale_x_continuous(labels = percent) + 
  xlab("% of 2017 Conservative promises kept\n(posterior density)") + 
  ylab("") + theme_minimal() + guides(colour = F) +
  theme(panel.grid.major.y = element_line(size=0.4, 
                                          color = "gray", linetype = 3),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.ticks.x = element_line(size=0.1, color = "black")) +
  geom_segment(aes(x=unweighted, xend = mean(props.total$Complete), y = 2.5, yend = 2.5),
               arrow = arrow(length = unit(0.3, "cm")), size = 1, colour = "#396a82") +
  # geom_point(aes(x = unweighted, y = 1), colour = "black", alpha = 0.1)  +
  geom_text(x = unweighted - 0.025, y = 3.2, 
            label = paste0("Unweighted\n", round(unweighted*100), "%")) + 
  geom_text(x = mean(props.total$Complete) - 0.02, y = 3.2, 
            label = paste0("Weighted\n", round(mean(props.total$Complete)*100), "%")) + 
  scale_color_gradient(low = "#dd250d", high = "#18a306") + 
  scale_y_discrete(drop = F)  + 
  geom_vline(xintercept = unweighted, linetype = 2) +
  geom_vline(xintercept = mean(props.total$Complete), linetype = 2)
density.all

saveForPub(density.all, file = "figures/density.all", height = 4, width = 7)


props.sq <- dtf(Complete = rowSums(propDone(imp.props, sq = c("yes"), pf= pf)))

unweighted.sq <- prop.table(table(pf$Done[pf$status_quo=="yes"]))["1"]
density.sq <- ggplot(data = props.sq, aes(x = Complete)) + 
  # scale_fill_manual(values = c("#3aa1fc", "#f49241")) +
  geom_density(colour = "#568aa3", alpha = 0.4, fill = "#3aa1fc") + 
  scale_x_continuous(labels = perc) + 
  xlab("% of 2017 Conservative promises kept\n(posterior density)") + 
  ylab("") + theme_minimal() + guides(colour = F) +
  theme(panel.grid.major.y = element_line(size=0.4, 
                                          color = "gray", linetype = 3),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.ticks.x = element_line(size=0.1, color = "black")) +
  geom_segment(aes(x=unweighted.sq, xend = mean(props.sq$Complete), y = 11.3, yend = 11.3),
               arrow = arrow(length = unit(0.3, "cm")), size = 1, colour = "#396a82") +
  # geom_point(aes(x = unweighted, y = 1), colour = "black", alpha = 0.1)  +
  geom_text(x = unweighted.sq - 0.007, y = 14.5, 
            label = paste0("Unweighted\n", round(unweighted.sq*100), "%")) + 
  geom_text(x = mean(props.sq$Complete) - 0.007, y = 14.5, 
            label = paste0("Weighted\n", round(mean(props.sq$Complete)*100), "%")) + 
  scale_color_gradient(low = "#dd250d", high = "#18a306") + 
  scale_y_discrete(drop = F, expand = c(0,0))  + 
  geom_vline(xintercept = unweighted.sq, linetype = 2) +
  geom_vline(xintercept = mean(props.sq$Complete), linetype = 2)+ geom_hline(size = 2, yintercept = 0)
density.sq

props.change <- dtf(Complete = rowSums(propDone(imp.props, sq = c("no"), pf= pf)))

unweighted.change <- prop.table(table(pf$Done[pf$status_quo=="no"]))["1"]

density.change <- ggplot(data = props.change, aes(x = Complete)) + 
  # scale_fill_manual(values = c("#3aa1fc", "#f49241")) +
  geom_density(colour = "#568aa3", alpha = 0.4, fill = "#3aa1fc") + 
  scale_x_continuous(labels = perc) + 
  xlab("% of 2017 Conservative promises kept\n(posterior density)") + 
  ylab("") + theme_minimal() + guides(colour = F) +
  theme(panel.grid.major.y = element_line(size=0.4, 
                                          color = "gray", linetype = 3),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.ticks.x = element_line(size=0.1, color = "black")) +
  geom_segment(aes(x=unweighted.change, xend = mean(props.change$Complete), y = 2.5, yend = 2.5),
               arrow = arrow(length = unit(0.3, "cm")), size = 1, colour = "#396a82") +
  # geom_point(aes(x = unweighted, y = 1), colour = "black", alpha = 0.1)  +
  geom_text(x = unweighted.change - 0.03, y = 3.2, 
            label = paste0("Unweighted\n", round(unweighted.change*100), "%")) + 
  geom_text(x = mean(props.change$Complete) - 0.03, y = 3.2,
            label = paste0("Weighted\n", round(mean(props.change$Complete)*100), "%")) + 
  scale_color_gradient(low = "#dd250d", high = "#18a306") + 
  scale_y_discrete(drop = F, expand = c(0,0))  + 
  geom_vline(xintercept = unweighted.change, linetype = 2) +
  geom_vline(xintercept = mean(props.change$Complete), linetype = 2) + geom_hline(size = 2, yintercept = 0)
density.change

posterior.densities <- arrangeGrob(density.change + ggtitle(paste0("Change promises (n=", sum(pf$status_quo=="no"), ")")), 
                                   density.sq + ggtitle(paste0("Status quo promises (n=", sum(pf$status_quo=="yes"), ")")), ncol = 2)
plot(posterior.densities)


getwd()

promise.complete.plots <- plotChangesFromImpProps(imp.props, imp.props.c, pf = pf)

done.draws <- dtf(Done_all_change = rowSums(propDone(imp.props = imp.props, "no", pf = pf)), 
                  Done_con_change = rowSums(propDone(imp.props = imp.props.c, "no", pf = pf)), 
                  Done_all_sq = rowSums(propDone(imp.props = imp.props, "yes", pf = pf)), 
                  Done_con_sq = rowSums(propDone(imp.props = imp.props.c, "yes", pf = pf)))



prop.table(table(done.draws$Done_con_change<prop.table(table(pf$statuscomb[pf$status_quo=="no"]))["Fulfilled"]))

density.all.done <- promise.complete.plots$change
density.all.sq<- promise.complete.plots$sq

library(mellonMisc)
library(grid)

hist(density.all.done$data$value[density.all.done$data$variable=="All"])


pf$Weight <- colMeans(imp.props)

density.ch.sq <- grid_arrange_shared_legend(density.all.done + ggtitle(paste0("Change promises (n=", table(pf$status_quo[!is.na(pf$statuscomb)])["no"], ")")),
                                            density.all.sq + ggtitle(paste0("Status quo promises (n=", table(pf$status_quo[!is.na(pf$statuscomb)])["yes"], ")")), ncol = 2)

plot(density.ch.sq)


promise.complete.plots.noeu <- plotChangesFromImpProps(imp.props = imp.props[, which(pf$MII_cat1!="europe")],
                                                       imp.props.c = imp.props.c[, which(pf$MII_cat1!="europe")],
                                                       pf = dtf(pf[which(pf$MII_cat1!="europe"), ]))

density.ch.sq.noeu <- grid_arrange_shared_legend(promise.complete.plots.noeu$change +
                                                   ggtitle(paste0("Change promises\n(n=", table(pf$status_quo[pf$MII_cat1!="europe"])["no"], ")")),
                                                 promise.complete.plots.noeu$sq +
                                                   ggtitle(paste0("Status quo promises\n(n=", table(pf$status_quo[pf$MII_cat1!="europe"])["yes"], ")")),
                                                 ncol = 2)

plot(density.ch.sq.noeu)


load(file = "data/intermediate/imp.simple.means.rda")


distributions.re.fmm <- dtf(RE = imp.simple.means, RE_FMM = colMeans(imp.samples))

distributions.re.fmm$FE <- fe$importFixedAll * 4

limits.dist = c(min(c(distributions.re.fmm$RE, distributions.re.fmm$RE_FMM, distributions.re.fmm$FE)), 
                max(c(distributions.re.fmm$RE, distributions.re.fmm$RE_FMM,distributions.re.fmm$FE)))


plot(distributions.re.fmm$FE, distributions.re.fmm$RE_FMM)
fe.hist <- ggplot(distributions.re.fmm, aes(x = FE)) +
  geom_histogram(colour = "black", fill = "lightblue") + 
  theme_bes()  + scale_y_continuous(expand = c(0,0)) + 
  scale_x_continuous(limits = limits.dist)+ xlab("Promise effect")

re.hist <- ggplot(distributions.re.fmm, aes(x = RE)) + geom_histogram(colour = "black", fill = "lightblue") + 
  theme_bes()  + scale_y_continuous(expand = c(0,0)) + 
  scale_x_continuous(limits = limits.dist)+ xlab("Promise effect")

re.fmm.hist <- ggplot(distributions.re.fmm, aes(x = RE_FMM)) + geom_histogram(colour = "black", fill = "lightblue") + 
  theme_bes()  + scale_y_continuous(expand = c(0,0)) + 
  scale_x_continuous(limits = limits.dist) + xlab("Promise effect")


comb.hist.re <- arrangeGrob(fe.hist + ggtitle("Fixed effects"), 
                            re.hist + ggtitle("Random effects"),
                            re.fmm.hist + ggtitle("Random effects/FMM"), ncol = 3)
plot(comb.hist.re)

saveForPub(comb.hist.re, file = "figures/comb.hist.re", width = 11, height = 4)


cor.c <- rep(NA, nrow(app.samples))

for(ii in 1:nrow(app.samples)) {
  if((ii / 100)==round(ii/100) ) {
    print(ii)
  }
  cor.c[ii] <- cor(as.matrix(dtf(imp = unlist(imp.samples.c[ii, ]), 
                                 app = unlist(app.samples.c[ii, ]))))[1,2]
}

cor <- rep(NA, nrow(app.samples))

for(ii in 1:nrow(app.samples)) {
  if((ii / 100)==round(ii/100) ) {
    print(ii)
  }
  cor[ii] <- cor(as.matrix(dtf(imp = unlist(imp.samples[ii, ]), 
                               app = unlist(app.samples[ii, ]))))[1,2]
}


library(matrixStats)


promise.summary.imp <- dtf(promise.index, 
                           Estimate = colMeans(imp.samples), 
                           est.median = colMedians(as.matrix(imp.samples)),
                           se = colSds(as.matrix(imp.samples)),
                           lci = colQuantiles(as.matrix(imp.samples), prob = 0.025), 
                           uci = colQuantiles(as.matrix(imp.samples), prob = 0.975), 
                           weight = colMedians(imp.props))

promise.summary.app <- dtf(promise.index, Estimate = colMeans(app.samples), 
                           est.median = colMedians(as.matrix(app.samples)),
                           se = colSds(as.matrix(app.samples)),
                           lci = colQuantiles(as.matrix(app.samples), prob = 0.025), 
                           uci = colQuantiles(as.matrix(app.samples), prob = 0.975))


promise.summary.imp.c <- dtf(promise.index, Estimate = colMeans(imp.samples.c), 
                             se = colSds(as.matrix(imp.samples)),
                             lci = colQuantiles(as.matrix(imp.samples.c), prob = 0.025), 
                             uci = colQuantiles(as.matrix(imp.samples.c), prob = 0.975),
                             weight = colMedians(imp.props.c))

promise.summary.app.c <- dtf(promise.index, Estimate = colMeans(app.samples.c), 
                             se = colSds(as.matrix(app.samples.c)),
                             lci = colQuantiles(as.matrix(app.samples.c), prob = 0.025), 
                             uci = colQuantiles(as.matrix(app.samples.c), prob = 0.975))



promise.summary.imp$weight <- prop.table(promise.summary.imp$weight)
promise.summary.imp.c$weight <- prop.table(promise.summary.imp.c$weight)
promise.summary.imp.c$est.z <- promise.summary.imp.c$Estimate / sd(promise.summary.imp.c$Estimate)
promise.summary.imp.c$weight <- prop.table(exp(promise.summary.imp.c$est.z))

any(is.na(match(promise.summary.imp.c$text, pf$Snippets)))

pf$Weight.c <- promise.summary.imp.c$weight[match(pf$Snippets, promise.summary.imp.c$text)]

pf$Imp.est <- promise.summary.imp$Estimate[match(pf$Snippets, promise.summary.imp$text)]
pf$Imp.lci <- promise.summary.imp$lci[match(pf$Snippets, promise.summary.imp$text)]
pf$Imp.uci <- promise.summary.imp$uci[match(pf$Snippets, promise.summary.imp$text)]


pf$Imp.est.c <- promise.summary.imp.c$Estimate[match(pf$Snippets, promise.summary.imp.c$text)]
pf$Imp.lci.c <- promise.summary.imp.c$lci[match(pf$Snippets, promise.summary.imp.c$text)]
pf$Imp.uci.c <- promise.summary.imp.c$uci[match(pf$Snippets, promise.summary.imp.c$text)]

pf$App.est <- promise.summary.app$Estimate[match(pf$Snippets, promise.summary.app$text)]
pf$App.lci <- promise.summary.app$lci[match(pf$Snippets, promise.summary.app$text)]
pf$App.uci <- promise.summary.app$uci[match(pf$Snippets, promise.summary.app$text)]


pf$App.est.c <- promise.summary.app.c$Estimate[match(pf$Snippets, promise.summary.app.c$text)]
pf$App.lci.c <- promise.summary.app.c$lci[match(pf$Snippets, promise.summary.app.c$text)]
pf$App.uci.c <- promise.summary.app.c$uci[match(pf$Snippets, promise.summary.app.c$text)]


pf$Imp.rank <- rank(pf$Imp.est)
pf$Imp.rank.c <- rank(pf$Imp.est.c)

pf$Imp.rank <- nrow(pf) - pf$Imp.rank + 1
pf$Imp.rank.c <- nrow(pf) - pf$Imp.rank.c + 1

pf$App.rank <- rank(pf$App.est)
pf$App.rank.c <- rank(pf$App.est.c)

pf$App.rank <- nrow(pf) - pf$App.rank + 1
pf$App.rank.c <- nrow(pf) - pf$App.rank.c + 1


pf$Imp.z <- pf$Imp.est/sd(pf$Imp.est)
pf$App.z <- pf$App.est/sd(pf$App.est)


pf$Imp.z.c <- pf$Imp.est.c/sd(pf$Imp.est.c)
pf$App.z.c <- pf$App.est.c/sd(pf$App.est.c)


est.plot <- ggplot(data= pf, aes(x = Imp.est)) + 
  geom_histogram(fill = "gray", colour = "black") + theme_bes() + 
  xlab("Random intercept")+ ylab("Count of promises") + 
  scale_y_continuous(expand = c(0,0)) + 
  ggtitle("Random Intercepts")

weight.plot <- ggplot(data= pf, aes(x = Weight)) + 
  geom_histogram(fill = "gray", colour = "black") + theme_bes() + 
  xlab("Exponentiated weight")+ ylab("Count of promises")+ scale_y_continuous(expand = c(0,0)) + 
  ggtitle("Exponentiated weight")

hist.comparisons <- arrangeGrob(est.plot, weight.plot, ncol = 2)
plot(hist.comparisons)

cor.mat <- cor(as.matrix(dtf(
  imp = promise.summary.imp$Estimate, 
  c.imp = promise.summary.imp.c$Estimate, 
  app = promise.summary.app$Estimate,
  c.app = promise.summary.app.c$Estimate)))

cor.mat[upper.tri(cor.mat)] <- NA
cor.mat["c.app", "c.imp"] <- median(cor.c)
cor.mat["app", "imp"] <- median(cor)

all.done <- sum(promise.summary.imp$weight[pf$Done==1 & pf$status_quo=="no"], na.rm = T)
all.not.done <- sum(promise.summary.imp$weight[!pf$Done==1 & pf$status_quo=="no"], na.rm = T)
prop.table(table(pf$Done))

cppg.done <- sum(promise.summary.imp$weight[pf$cppg & pf$Done==1& pf$status_quo=="no"], na.rm = T)
cppg.not.done <- sum(promise.summary.imp$weight[pf$cppg & !pf$Done==1& pf$status_quo=="no"], na.rm = T)

save(promise.summary.imp.c, file="data/intermediate/promise.summary.imp.c.rda")
save(promise.summary.imp, file="data/intermediate/promise.summary.imp.rda")
save(promise.summary.app.c, file="data/intermediate/promise.summary.app.c.rda")
save(promise.summary.app, file="data/intermediate/promise.summary.app.rda")





save(pf, file = "data/intermediate/pf_results.rda")



#### OUTPUT ####

saveForPub(density.all.done, file = "figures/density.all.done", width = 7, height = 5)


sort(tapply(fe$importWtAll, fe$MII_cat1, sum))
sort(tapply(fe$importWtCon, fe$MII_cat1, sum))


# Table 1 Central promises: combined coefficient centered at zero
# (centrality and approval), rank (centrality and approval) and weight 
# (centrality only). 
pf %>% slct(Promise = Snippets, Issue = agreed_cat_string, 
            Status = statuscomb, 
            Imp.est, Imp.rank, Imp.weight = Weight, 
            App.est, App.rank) %>% arrange(Imp.rank) %>% 
  mutate(Imp.weight = round(Imp.weight * 100,1), 
         App.est = round(App.est, 1), 
         Imp.est = round(Imp.est, 1), 
         Status = Hmisc::capitalize(Status), 
         Issue =  Hmisc::capitalize(Issue)) %>% 
  mutate(Status = recode(Status, "Not fulfilled" = "Broken"))  %>% 
  filter(Imp.rank %in% 1:20) %>%
  write.csv(file = "tables/top20promises.csv")


# First, we examine the correlation between approval and centrality, and find they are uncorrelated (r=0.08).  
cor.mat["app", "imp"]

# The most central promise (“Leave the EU”) receives 9% of all the weight. 
# ...
# The second largest promise to “Reduce annual net migration to under 100,000”
# receives 8% weight. 
pf %>% arrange(desc(Weight)) %>% head(5) %>% slct(Snippets, Weight)

# EU-related promises make up 35% of weighted promises
head(sort(tapply(pf$Weight, pf$agreed_cat_string, sum), decreasing = T),10)
# , compared with just 5% if we count all promises equally.
# (intro) If we take the standard approach of treating each promise equally,
# just 5% of Conservative promises related to Brexit
prop.table(table(pf$agreed_cat_string))["europe"]

# In particular, the combined “other” category (all 24 remaining issue areas) 
# which contains low weighted promises such as “Use 3rd party identification 
# to make accessing online government services more secure” and “Sanction firms
# that don’t obey government internet regulations” falls dramatically from 61% 
# of the Conservatives’ agenda when equally weighted to 27% using the centrality weights. 

prop.table(table(pf$cat_small))["other"]
sort(tapply(pf$Weight, pf$cat_small, sum), decreasing = T)

# Abstract: Weighting increases the centrality of EU promises sevenfold, immigration 
# promises sixfold, and reduces the centrality of miscellaneous administrative 
# promises by more than half. 
# (intro) , Brexit’s importance increases sevenfold.
tapply(pf$Weight, pf$cat_small, sum) / prop.table(table(pf$cat_small))


# Figure 4. Posterior distributions of the percentage of Conservative manifesto 
# in issue areas. Arrows indicate the difference between equally-weighted
# proportion of promises and the median weight falling into that category. 
# Promise categories separately ordered by median value. Other category includes
# combined weight of all promises not included in 7 main categories listed.
saveForPub(joy.all.iss, file = "figures/joy.all.iss", height = 4, width = 7)

# We find the equally-weighted pledge fulfilment of the government was 69%,
# a rate comparable to other pledge fulfilment studies 
prop.table(table(pf$statuscomb))

# However, when we weight the promise fulfilment according to the conjoint 
# experiment, we find only 48% of promises were fulfilled. 
prop.table(tapply(pf$Weight, pf$statuscomb, sum))

# abstract: Centrality weighting reduces our assessment of Conservative 
#...
# The Conservatives’ pledge fulfilment rate falls by 21 percentage points after 
# weighting, equivalent to a two standard deviation fall in fulfilment
# promise-keeping by 21 percentage points (1.3 standard deviations of typical
# promise completion rates found in comparative studies). 
# ...
# The reduction of pledge fulfilment by 21 percentage points 
# ...
# (conclusions) Accounting for pledge centrality when assessing the 2017 
# Conservative party’s manifesto reduces our assessment of its promise
# completion by 21 percentage points. 
# is 131% of the size
diffs <- prop.table(tapply(pf$Weight, pf$statuscomb, sum)) - prop.table(table(pf$statuscomb)) 
diffs["Fulfilled"]
load(file = "data/intermediate/sd.cppg.fulfil.rda")
diffs["Fulfilled"] / sd.cppg.fulfil


# Footnote: FE models give 46% fulfilment for change promises.
# Footnote: FE models give 95% fulfilment for status quo promises.
prop.table(tapply(fe$importWtAll, list(sq = fe$status_quo, status = fe$statuscomb), sum), 1)

#   FE models give EU promises a 24% weight. 
sort(tapply(fe$importWtAll, fe$agreed_cat_string, sum))

# Figure 5 separates completion by whether a promise was to change (216 promises) 
# or maintain the status quo (39 promises). 
table(pf$status_quo[!is.na(pf$statuscomb)])


# For promises to change the status quo, we see a substantial difference when we 
# use pledge-centrality weights. Figure 5 shows that the Conservatives fulfilled 
# 40% of their promises to change the status quo when we use these weights, 
# compared with 64% when weighting all promises equally.  
# ...
# The pledge-centrality weights suggest the Conservatives fulfilled 96% of their
# promises to maintain the status quo, whereas this figure was only 92% when 
# weighting promises equally
prop.table(table(statuscomb = pf$statuscomb, status_quo = pf$status_quo), 2)
prop.table(tapply(pf$Weight,list(statuscomb = pf$statuscomb, status_quo = pf$status_quo), sum), 2)


# Our results show a >99.9% probability that pledge-centrality weighting lowers 
# the proportion of change promises completed. 
prop.table(table(done.draws$Done_all_change<prop.table(table(pf$statuscomb[pf$status_quo=="no"]))["Fulfilled"]))


# Figure 5.
# Posterior distributions of Conservative promises fulfilled (%) using equal 
# weights and weights derived RE/FMM models of the conjoint experiment. 
# Arrows indicate how our assessment of fulfillment moves from equal weighting 
# to the median of the posterior distribution of weighted completion.
saveForPub(posterior.densities, file = "figures/posterior.densities", width = 12, height = 5)

prop.table(table(statuscomb = pf$statuscomb, status_quo = pf$status_quo), 2) - 
  prop.table(tapply(pf$Weight,list(statuscomb = pf$statuscomb, status_quo = pf$status_quo), sum), 2)

# If May had passed her Brexit bill through parliament, the government’s 
# completion rate would have increased by 23 percentage points. 
sum(pf$Weight[pf$Snippets %in% c("Leave the EU",  "Leave the EU's Single Market",
                                 "Leave the EU's Customs Union" )])


#### Appendix A ####
# In all, we removed 37 promises. 
table(pf$cppg)

# The tighter definition slightly decreases the proportion of promises to change
# the status quo completed when using weights (38%) compared to weighting using 
# the wider definition (40%). 
prop.table(c(all.done, all.not.done))
prop.table(c(cppg.done, cppg.not.done))


# Using the tighter definition, the median percentage of promise weight relating
# to the EU is slightly higher (37%) 
prop.table(tapply(pf$Weight[pf$cppg], pf$cat_small[pf$cppg], sum))
# than when using the full definition (35%). 
prop.table(tapply(pf$Weight, pf$cat_small, sum))



#### Appendix E ####

# Appendix E: Table 1 Centrality and approval coefficients, rank, 
# and weight for all respondents and 2017 Conservative voters
pf %>% slct(Promise = Snippets, Issue = agreed_cat_string, 
            Status = statuscomb, 
            Imp.est, Imp.est.c, 
            Imp.rank, Imp.rank.c, 
            Imp.weight = Weight, Imp.weight.c = Weight.c, 
            App.est, App.est.c, 
            App.rank, App.rank.c) %>% 
  arrange(Imp.rank) %>% 
  mutate(Imp.weight = round(Imp.weight * 100,1), 
         Imp.weight.c = round(Imp.weight.c * 100,1), 
         App.est = round(App.est, 1), 
         App.est.c = round(App.est.c, 1), 
         Imp.est = round(Imp.est, 1), 
         Imp.est.c = round(Imp.est.c, 1), 
         Status = Hmisc::capitalize(Status), 
         Issue =  Hmisc::capitalize(Issue)) %>% 
  mutate(Status = recode(Status, "Not fulfilled" = "Broken"))  %>%
  write.csv(file = "tables/allpromises.csv") 


# Figure 8 Histograms of the random intercepts (left panel) and weights
# (right panel) of the 257 promises after using our modelling approach.
saveForPub(hist.comparisons, file = "figures/hist.comparisons", height = 4, width = 10)



#### Appendix F ####

# Both models agree that EU-related promises make up the most important issue 
# area, making up 35% of weighted promises for the all voters model and
prop.table(tapply(pf$Weight, pf$cat_small, sum))

# 41% for Conservative voters.
prop.table(tapply(pf$Weight.c, pf$cat_small, sum))





# 61% of the Conservatives’ agenda when equally weighted to 27% using the 
# centrality weights from the all voters model, and only 18% when
# weighted using Conservative voters’ responses.
prop.table(table(pf$cat_small))["other"]
prop.table(tapply(pf$Weight, pf$cat_small, sum))["other"]
prop.table(tapply(pf$Weight.c, pf$cat_small, sum))["other"]

# Using the all voters model we find only 48% of promises were fulfilled, 
# using the Conservatives only model this number is lowered further, 
# with a fulfilment rate of only 44%.
prop.table(tapply(pf$Weight, pf$statuscomb, sum))
prop.table(tapply(pf$Weight.c, pf$statuscomb, sum))


# Using the all-voters weights, we find that the Conservatives partially or 
# fully fulfilled 40% of their promises to change the status quo. 
prop.table(tapply(pf$Weight[!is.na(pf$statuscomb)], 
                  list(pf$status_quo[!is.na(pf$statuscomb)],
                       pf$statuscomb[!is.na(pf$statuscomb)]), sum), 1)
# Using the Conservatives-only weights we find a completion rate of 38%. 
prop.table(tapply(pf$Weight.c[!is.na(pf$statuscomb)], 
                  list(pf$status_quo[!is.na(pf$statuscomb)],
                       pf$statuscomb[!is.na(pf$statuscomb)]), sum), 1)

# For promises to keep the status quo we find a fulfilment rate of 96% using
# the all-voters model
prop.table(tapply(pf$Weight[!is.na(pf$statuscomb)], 
                  list(pf$status_quo[!is.na(pf$statuscomb)],
                       pf$statuscomb[!is.na(pf$statuscomb)]), sum), 1)

# and 96% for the Conservatives-only weights. 
prop.table(tapply(pf$Weight.c[!is.na(pf$statuscomb)], 
                  list(pf$status_quo[!is.na(pf$statuscomb)],
                       pf$statuscomb[!is.na(pf$statuscomb)]), sum), 1)

# nine out of ten most central promises are the same promises in both models
table(pf$Imp.rank %in% 1:10, pf$Imp.rank.c %in% 1:10)


# Figure 9 Posterior distributions of the percentage of Conservative manifesto 
# in issue areas for all respondents and Conservatives voters only. Arrows
# indicate the difference between unweighted proportion of promises and the
# median weight falling into that category. Promise categories separately 
# ordered by median value. Other category includes combined weight of all
# promises not included in 7 main categories listed.
saveForPub(joy2, file = "figures/joy2", height = 9, width = 7)

# Figure 10 Posterior distributions of percentage of Conservative party promises
# fulfilled using equal weights and weights derived RE/FMM models of the
# conjoint experiment for all respondents and Conservatives voters only. 
# Arrows indicate the movement of our assessment of fulfillment from equal 
# weighting to the median of the posterior distribution of weighted completion.

saveForPub(density.ch.sq, file = "figures/density.ch.sq", width = 12, height = 4)

#### Appendix I: ####

# promise weighting still substantially reduces the percentage of the
# Conservatives’ agenda to change the status quo that has been completed 
# (18 percentage point reduction for all respondents and a 33 percentage point 
# reduction for 2017 Conservative voters). 
promise.complete.plots.noeu$change$data$Unweighted[1] - 
  mean(promise.complete.plots.noeu$change$data$value[promise.complete.plots.noeu$change$data$variable=="All"]) 

promise.complete.plots.noeu$change$data$Unweighted[1] - 
  mean(promise.complete.plots.noeu$change$data$value[promise.complete.plots.noeu$change$data$variable=="2017 Conservatives"]) 





# Appendix I: Figure 12 :  Posterior distributions of the percentage of Conservative
# manifesto in issue areas according to all respondents and 2017 Conservative 
# voters excluding Europe promises. Arrows indicate the difference between 
# unweighted proportion of promises and the median weight falling into that category. 
# Promise categories separately ordered by median value for all respondents 
# and 2017 Conservatives. Other category includes combined weight of all
# promises not included in 7 main categories listed.
saveForPub(joy.all.iss.noeu, file = "figures/joy.all.issnoeu", height = 4, width = 7)

# Appendix I: Figure 13 Posterior distributions of percentage of Conservative 
# party promises fulfilled using equal weights and weights derived RE/FMM models 
# of the conjoint experiment excluding Europe promises. Arrows indicate the 
# movement of our assessment of fulfillment from equal weighting to the median
# of the posterior distribution of weighted completion.
saveForPub(density.ch.sq.noeu, file = "figures/density.ch.sq.noeu", width = 12, height = 4)


# Appendix F: Table 2 Correlation coefficients between random intercepts of 
# approval and centrality for full sample/2017 Conservative sample (n=257).
# Correlations within a model are median correlation across MCMC iterations.

round(cor.mat,2)
write.csv(cor.mat, file = "tables/cor_mat.csv")




